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ABSTRACT 

In this paper we present an algorithm to address the pre- 
decessor problem of feed-forward Boolean networks. We 
propose an probabilistic algorithm, which solves this prob- 
lem in linear time with respect to the number of nodes in 
the network. Finally, we evaluate our algorithm for ran- 
dom Boolean networks and the regulatory network of Es- 
cherichia coli. 

1. INTRODUCTION 

In systems and computational biology Boolean networks 
(BN) are widely used to model regulative dependencies of 
organisms [T T]. We consider networks, which map a set 
of environmental conditions to the presence of proteins 
and finally to actual chemical reactions, which are often 
modeled as fluxes of flux-balance analysis pl|. Hence, 
these networks are used to make in silico predictions of 
behavior of organisms in a certain environment H. 

In this paper we address the inverse problem, i.e., we 
want to predict environmental conditions that allow cer- 
tain reactions to take place, and others not. Hence, in 
general, we need to find a set of possible inputs that lead 
to a given output. This so called predecessor problem 
or preimage problem has been addressed by Wuensche in 
Q and has been shown to NP-hard in general f6l, which 
makes it infeasible to solve it for large networks. In [TJ an 
algorithm with reduced complexity for BNs with canaliz- 
ing Boolean functions has been introduced. However, the 
problem is infeasible under certain conditions. Both al- 
gorithms are designed to find the whole set of preimages, 
i.e., all inputs to the BN with lead to a certain, desired, 
output. 

In some applications, knowledge of the whole preim- 
age set is not important, merely it can be sufficient to know 
a subset of the preimage set. Here, we propose an proba- 
bilistic algorithm, which solves this problem in linear time 
with respect to the number of nodes in the network, based 
on a variation of the well known Sum-Product- Algorithm 
m, which is used for a variety of tasks, including de- 
coding error correction codes in communication engineer- 
ingL9J. 

2. BOOLEAN NETWORKS AND MAIN IDEA 

We consider networks like shown in Figure [1] mapping 
the values of the N in-nodes I = {1, 2, 3} to the M out- 




Figure 1 . Example of a Feed-Forward Network 

nodes = {12, 13, 14, 15, 16}, i.e., we can represent this 
BN as a function mapping the N input values uniquely to 
the M output values: 



f :{0,ir ^{0,1} 



M 



The network itself consists of n nodes, and a set of di- 
rected edges connecting these nodes. Each node i has a 
certain state, which can be either zero or one, represented 
by a variable Xi. Its value is determined by evaluating 
a Boolean function (BF) fi. Further, lets define the set 
n{fj) as the incoming nodes of node j. For example in 
Figure [T] ^(/s) = {1,3}. The BF fj is a function map- 
ping kj = \n{fj)\ values of {0, 1}'^ to {0, 1}, where k is 
also called the in-degree of node j. The number of edges 
emerging from a node is called out-degree. 

Given a vector of input values x G {0,1}^, x = 
{xi, X2, ■ ■ ■ ,xn) the corresponding output of f is y = 
f(x),y G {0,1}*^. In general there does not exist a 
unique inverse function f Instead the cardinality of the 
set := {x : f (x) = y} will be larger one. We call fly 
the set of preimages of y. 

In this paper we are interested to find at least parts 
of fly. Suppose there is a probability distribution Py on 
{0,1}^ such that 



Py{x} = \y 







if X G 
else 



If we knew the probability distribution Py, we would have 
solved the problem. But as explained, this is too difficult 
in general. Our main idea now is to approximate Py by the 
product of the marginal distributions Pi on the individual 



Xi, i.e., 

N 

1=1 

as the well-known Sum-Product algorithm can be used to 
compute the marginals efficiently. If the approximation is 
good enough sampling out the product of the marginals 
will yield an element in ily with reasonable probability. 

3. PROPOSED ALGORITHM 

In this section we will first discuss the basic principles 
of factor graphs and the Sum-Product Algorithm (Section 
13.1b . Then we will describe the BN as factor graph in 
Section [J!2] and will formulate the actual algorithm to find 
the marginals in Section [33] Finally, the sampling is de- 
scribed. 

3.1. Factor Graphs and Sum-Product Algorithm 

Assume some function 5 (xi, . . . , a;„) defined on some do- 
main A", which can be factorized in m local functions 
hj,j e H := {1,2,...,™}, i.e., 

g{xi, . . . ,x„) = 

j 

where Xj is the subset of [n] containing the argument 
of hj. We can then define a factor graph |l8] as a bi- 
partite graph consisting of n nodes representing variables 
{a^i, . . . ,Xn} (variable nodes) and of m nodes represent- 
ing functions {/i, . . . fm} (function node). Edges only 
exist between a function node and a variable node if and 
only if Xi is an input to function fj . 

The marginal function gi{xi) is defined as lH) 

g^{x^) = ^ g{xi,...,Xn), 
~{xi} 

where Y.^{xi} 9{xi, . . . , a;„) is defined as 
^ g{xi,...,Xn) 

In general the computation of the is difficult, but due 
to the factorization of g the task can be efficiently solved 
using the the so called Sum-Product algorithm [8|. The 
algorithm iteratively passes messages between the nodes 
of the graph. At each iteration the messages /x are sent 
from the function nodes to the variable nodes, containing 
the corresponding marginal function of the local function. 
These messages are computed as follows JS): 
function to variable node: 

-{an} \ y(in(h)\{x} ) 

where n{i) give the set of neighboring nodes of node i. 



At the variable nodes, these messages are then com- 
bined to a marginal function A and sent back to the func- 
tion nodes ISj: 

variable to function node: 

\x^h{x) = W pLq^x{x). 

qen(x)\{h} 

3.2. The Boolean Network as Factor Graph. 

We apply the concept of factor graphs to BNs. Each node 
in the network represents one variable Xi G {0, 1}, i e [n] 
of the factor graph, hence we have n variable nodes. Each 
BE fj of the BN {j G [n] \ I) is a function node and is 
connected to the node j and the incoming nodes n{fj). 
Lets to define Xj as the variables of the incoming nodes 
of node j, i.e. the argument of the BN fj. Eurther, we 
define X - ' as X, without the node i. 

Einally, if we consider the variables as each node as 
random variables, we have a common distribution of all 
variables nodes described by the density function, 

gX \ ,...,Xn {Xl ,...5X7^) = {Xl , . . . , Xyj ) , 

Eor sake of readability we will omit the subscript of the 
density function, if they are obvious from context. We 
are interested in finding the marginal distributions of the 
in-nodes, which can be described by the density functions 

gxiixi) = ^gxi,...,x„ixi, ...,x„) Vi e L 

^Xi 

This problem is an instance of the problem described in in 
Section im hence, we apply the same methods here. 

3.2.1. Update Rule: function to variable node 

If we focus on one function node j e [j^] \ I there exists a 
common distribution of all variables relevant for this node. 
Namely, these relevant variables are the ones located in 
Xj of the BE fj, and the value of node j. We can write 
the density of this distribution as: 

p{xj,Xj). 

Lets define h{fj) as the set of indices of the input nodes 
of the BF fj. 

We need to send the local marginal distribution of each 
variablei G {j}Un(/j) back to the variable node, or more 
formal: 

lij^,{x,) ^ p{xj,Xj)= p(xj,x^,xf) 

~{a:i} ~{xi} 

(1) 

\li — i , i.e. if the message is designated for the node con- 
taining the output of the BE, the density of the marginal 
distribution becomes: 

~{xj} 



which is the probability distribution of the functions out- 
put. We can assume that the elements of Xj are pairwise 
independent, hence, we can write: 

where A; is the probability distribution of node j and is 
defined in Section [3.2.2l 

In the other cases, i.e., i ^ j, Eq. ([TJ becomes: 

~{xi} 

We still can assume that the elements of X^ are pairwise 
independent, hence, we can write: 

p{x,Mfj)\x^^p{^,\xf)-p{xf) 

^p{x,\xf) n \i{xi). 

l&Hf,}\{3} 

If the Boolean functions output xj = fj{Xj) is al- 
ready completely determined by Xj i.e., if the variable 
Xi has no influence on the output for this particular choice 
of the other variables, we assume Xi to be uniformly dis- 
tributed: 



p{x^\xj,X^j"'>) = ^ 
and since xj is completely determined by X 



1 



p{x,,xf)^ n ^'(^')- 

Otherwise, Xi is totally determined by Xj and the other 
variables, i.e., Xi is or 1 depending on BE Hence, we 
can write 

p{xi\xj,n{fj) \ xi) ^ p-^.{f{X^j'\x,) = Xj), 

where p^^ Xi) — Xj) is either or 1. Eurther we 

can assume Xj independent of ' , hence, 

p{x,,xf)^\,{xj) n ^'(^')- 

'ef.(/,)\{i} 
Einally, we can summarize for i ^ j: 

~{x,} ieft(/,)\b} 



with 



(0 



K^^j 



Xj {xj ) , else 
3.2.2. Update Rule: variable to function node 

The update rule is the same for all variable nodes j G [n] 
and is independent of the function node to which they are 
directed. 

^ji^s) ^W_l^i~^]ixj)^ (3) 

where Sj is the set of all function nodes, which have node j 
as input. 



3.3. Finding the Input Distributions 

In our algorithm, we use the well known log-likelihood 
ratio (LLR) to represent the probability distribution of bi- 
nary variables ifTOl . It is defined as: 



Lx = In 



p{x = 0) 
p{x = !)■ 



(4) 



A scheme of the algorithm is given in Algorithm[T] 

The probability distribution of each node j G [n] at it- 
eration t is given as L^p and are initialized with L^"'' = 0, 
which is equivalent to the uniform distribution. Then we 
set the LLRs for the out-nodes to either — oo or +oo de- 
pending on the desired output y of the BN. At each itera- 
tion the algorithm can be split in two steps. The first step 
iterates over all function nodes j G [n] \ I and all input 
variables i £ n{fj) calculating the LLR using Eq. 
(|2| and Eq. 

In the second step we update all variables-nodes, where 
the LLRs Lj represents the distributions Xj and, hence, 
the product of Eq. [3]becomes a summation. Please note, 
that the LLR of the previous iteration is also added to the 
sum, in order to prevent rapid changes of the distributions. 

After performing a certain number of iterations tmax, 
the desired marginal distributions of the input variables 
are found. 

Algorithm 1 



Initialize L 



(0) 



for all nodes 



Set the desired LLRs of the out-nodes, i.e., ij"^ is either 
— oo or +00, for all out-nodes j G 0. 
t=0 

repeat 

t=tH-l 

for each non-in-node j G [»^] \ I do 

for each input variable i G n{fj) do 

calculate Lj^i using Eq. (|2]i and Eq. (|4]l 

end for 
end for 

for each non-out-node v do 



end for 

until maximum number of iterations reached 



3.4. Sampling 

The sampling part of our approach is straight forward. Us- 
ing the marginal distributions t^^*™"^) ^ j G I we randomly 
draw vectors x and check if they fulfill y ~ /(x). If so, 
they are added to the set fly. This procedure is repeated 
for a certain number of samples. 

4. SIMULATION RESULTS AND DISCUSSION 

We tested our algorithm with randomly generated networks 
and the regulatory network of Escherichia coli (E-coli) 
||2l. The random networks consist of 2400 nodes with 
N = 200 and M 1200. We have chosen the BEs from: 




10 20 30 

imax 



Figure 2. Similarity of y and y verses tmax 

• all functions with A; < 15 (Type A) 

• unate, i.e. locally monotone, functions with fc < 15 
(TypeB) 

After generating a network we draw a certain number T 
of uniformly distributed input vectors x and obtain y = 
f (x). For each y we applied then Algorithm [T] to obtain 
the marginal distributions t^^*™"^) j g | 

To investigate the convergence behavior with respect 
to t„iax and we first apply hard-decision to evaluate a good 
choice for tmax, i-C-, we generate an estimate x by setting 

, _ Jo ifij*— ) > 
""'-{l ifLf— )<0 

Then we evaluate the network y — f (x), and measure 
the similarity between y and y by counting the equal en- 
tries and divide them by the length of y. We did so for 
100 networks of Type A and B, and set T = 100. The 
averaged results can be seen in Figure |2] One can see, that 
for tmax > 14 there is almost no improvement in the sim- 
ilarity. This number is equal to two times the number of 
nodes between input and output, i.e., it seems to be suf- 
ficient that the messages travel once through the network 
and back. Thus, the following simulations have been per- 
form setting tmax — 14. 

Next, we apply sampling as described in Section l?!4l 
We did so for 100 different networks of Type A and B, 
and the E-coli network. For each random network we did 
T = 100 runs, for E-coli T = 1000. The results can be 
viewed in Table [T] We depict the percentage of solved 
networks, i.e. the portion of networks we found at least 
one valid x G fiy. Further, we give the average number 
of valid x and the average number of unique x. 

One can see from the results, that in general for most 
networks and ys at least one preimage can be found. It is 
worth mentioning, that for the E-coli network every sam- 
pled solution was unique. This is due to the fact, that 
there exist a few inputs, who completely determine the 



network 


num of samples 


solved 


valid 


unique 


Type A 


1000 


89% 


608.81 


4.43 


TypeB 


1000 


95.9% 


270.74 


68.60 


E-coli 


1000 


98.6% 


193.3 


193.3 



Table 1 . Simulation results for different networks 



output. The other input variables have then no influence 
and hence a marginal distribution of 0.5. Further, the re- 
sults for the network of type B are much better than for 
type A. It seems that the marginal distributions for unate 
functions give better estimation of the actual distribution 
than the marginal distributions for non-unate functions. 
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